data{
  int<lower=0> N;
  int<lower=1> K; //number of category
  int<lower=1> M; // number of covariates
  matrix[N, M] covariates;
  int<lower=0, upper=1> treatment[N];
  int<lower=1,upper=K> Y[N];
}
parameters{
  vector[M] betaX; // beta for covariates
  real beta; // beta for treatment
  ordered[K-1] c;
}
model{
  for(i in 1:N){
    Y[i] ~ ordered_logistic(covariates[i,] * betaX + beta*treatment[i], c);
  }

}
generated quantities{
	real T[K];
	real C[K];
	real diff[K];
  real diff_D[2];  // dichotomized
  real diff_DwN[2];  // dichotomoized with neutral

	for (k in 1:K) {
		T[k]=0;
		C[k]=0;
		diff[k]=0;
	}

  for(k in 1:2) {
    diff_D[k] = 0;
    diff_DwN[k] = 0;
  }

	for(i in 1:N){

		T[1] = 1 - inv_logit(covariates[i,] * betaX + beta*1 - c[1]);
		C[1] = 1 - inv_logit(covariates[i,] * betaX + beta*0 - c[1]);	
		diff[1] = diff[1] + T[1] - C[1];

		for(s in 2:K-1){
			T[s] =  inv_logit(covariates[i,] * betaX + beta*1 - c[s-1]) - inv_logit(covariates[i,] * betaX + beta*1 - c[s]);
			C[s] =  inv_logit(covariates[i,] * betaX + beta*0 - c[s-1]) - inv_logit(covariates[i,] * betaX + beta*0 - c[s]);
			diff[s] = diff[s] + T[s] - C[s];
		}

		T[K] =  inv_logit(covariates[i,] * betaX + beta*1 - c[K-1]);
		C[K] =  inv_logit(covariates[i,] * betaX + beta*0 - c[K-1]);
		diff[K] = diff[K] + T[K] - C[K];

    diff_D[1] = diff[1] + diff[2];
    diff_D[2] = diff[4] + diff[5];

    diff_DwN[1] = diff[1] + diff[2];
    diff_DwN[2] = diff[3] + diff[4] + diff[5];

	}

	for(k in 1:K) {
		diff[k] = diff[k]/N;
  }

  for (k in 1:2) {
    diff_D[k] = diff_D[k]/N; 
    diff_DwN[k] = diff_DwN[k]/N; 
  }

}
